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Abstract 

The general number field sieve (GNFS) is the most efficient algorithm 
known for factoring large integers. It consists of several stages, the first 
one being polynomial selection. The quality of the chosen polynomials in 
polynomial selection can be modelled in terms of size and root properties. 
In this paper, we describe some algorithms for selecting polynomials with 
very good root properties. 

1 The general number field sieve 

The general number field sieve [TT] is the most efficient algorithm known 
for factoring large integers. It consists of several stages including polyno- 
mial selection, sieving, filtering, linear algebra and finding square roots. 

Let n be the integer to be factored. The number field sieve starts 
by choosing two irreducible and coprime polynomials f(x) and g(x) over 
Z which share a common root m modulo n. In practice, the notations 
F(x, y) and G(x, y) for the homogenized polynomials corresponding to / 
and g are often used. We want to find many coprime pairs (a, b) € 1? 
such that the polynomials values F(a, b) and G(a, b) are simultaneously 
smooth with respect to some upper bound B. An integer is smooth with 
respect to bound B (or B-smooth) if none of its prime factors are larger 
than B. The lattice sieving [16] and line sieving [4] are commonly used 
to identify such pairs (a,b). The running-time of sieving depends on the 
quality of the chosen polynomials in polynomial selection, hence many 
polynomial pairs will be generated and optimized in order to produce a 
best one. 
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This paper discusses algorithms for root optimization in polynomial 
selection in the number field sieve. We mainly focus on polynomial selec- 
tion with two polynomials, one of which is a linear polynomial. 

2 Polynomial selection 

For large integers, most methods [4j [9] 1101 1131 114] use a linear polynomial 
for g(x) and a quintic or sextic polynomial for f(x). Let f(x) = X/»=o CiX * 
and g(x) = m^x — mi. The standard method to generate such polynomial 
pairs is to expand n in base-(mi, 7712) so n = X^f=o c s m i m 2~ l ■ 

The running-time of sieving depends on the smoothness of the poly- 
nomial values \F(a, b)\ and |G(a,6)|. Let ^(x, x 1 ^) be the number of 
^/"-smooth integers below x for some it. The Dickman-de Bruijn func- 
tion p(u) [7] is often used to estimate ^(x,x x ^ u ). It can be shown that 

hm — - — p{u). 

a?— >oo X 

The Dickman-de Bruijn function satisfies the differential equation 

up'(u) + p{u - 1) = 0, p{u) = 1 for < u < 1. 

It may be shown that p satisfies the asymptotic estimate 

log(p(u)) = —(1 + o(l))u\ogu as u — > 00. 

For practical purposes, the frequency of smooth numbers can be approxi- 
mated by the Canfield-Erdos-Pomerance theorem, which can for example 
be stated as follows [8]. 

Theorem 2.1. For any fixed e > 0, we have 

V(x,x 1/U )=xu-^ 1+ ° W) 

as a; 1 /" and u tends to infinity, uniformly in the region x > u"'' 1_e . 

It is desirable that the polynomial pair can produce many smooth in- 
tegers across the sieve region. This heuristically requires that the size of 
polynomial values is small in general. In addition, one can choose an alge- 
braic polynomial f(x) which has many roots modulo small prime powers. 
Such a choice is driven by inheritance of practices which already date 
back to the CFRAC era, where suitable multipliers were chosen precisely 
in order to optimize this very property [121 117] . Then the polynomial 
values are likely to be divisible by small prime powers. This may in- 
crease the smoothness chance for polynomial values. We describe some 
methods [5][]3] to estimate and compare the quality of polynomials. 

2.1 Sieving test 

A sieving experiment over short intervals is a relatively accurate method to 
compare polynomial pairs. It is often used to compare several polynomial 
candidates in the final stage of the polynomial selection. Ekkelkamp [5] 
also described a method for predicting the number of relations needed 
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in the sieving. The method conducts a short sieving test and simulates 
relations based on the test results. Experiments show that the prediction 
of the number of relations is within 2% of the number of relations needed 
in the actual factorization. 

2.2 Size property 

Let (a, b) be pairs of relatively prime integers in the sieving region f2. For 
the moment, we assume that a rectangular sieving region is used where 
\a\ < U and < b < U. We also assume that polynomial values \F(a, 6)| 
and \G(a, b)\ behave like random integers of similar size. The number of 
sieving reports (coprime pairs that lead to smooth polynomial values) can 
be approximated by 

n 

The multiplier Q/ty 2 accounts for the probability of a, b being relatively 
prime. 

Since G is a linear polynomial, we may assume that log(| G(a, b)\) does 
not vary much across the sieving region. A simplified approximation to 
compare polynomials (ignoring the constant multiplier) is to compare 

The base- (mi, 7712) expansion [5J [TU] gives polynomials whose coeffi- 
cients are 0(n 1 ^ d+1 ' > ). The leading coefficients a and Cd-i are much 
smaller than 7i 1//(d+1) . The coefficient c d -2 is slightly smaller than n 1/(d+1) . 
For such polynomials, it is often better to use a skewed sieving region 
where the sieving bounds for a, b have ratio s, while keeping the area of 
the sieving region 2U 2 . The sieving bounds become \a\ < U^/s and < 
b < U l\fs. Each monomial in the polynomial is bounded by CiU d s l ~ d / 2 . 

In the integral |T]), computing p is time-consuming, especially if there 
are many candidates. We can use some coarser approximations. Since 
p(u) is a decreasing function of u, we want to choose a polynomial pair 
such that the size of |.F(a, 6)| and |G(o, 6)| is small on average over all 
(a, 6). This roughly requires that the coefficients of the polynomials are 
small in absolute value. We can compare polynomials by the logarithm of 
a L 2 -norm for polynomial F(x,y) by 

~log(s- d j J ^F 2 (xs,y)dxdy) . (2) 

where s is the skewness of sieving region. Polynomials which minimize 
the expression [5] are expected to be better than others. 

2.3 Root property 

If a polynomial f(x) has many roots modulo small prime powers, the 
polynomial values may behave more smoothly than random integers of 
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about the same size. Boender, Brent, Montgomery and Murphy [3] 1131 
1X41 1 1 5] described some quantitative measures of this effect (root property). 

Let p be a fixed prime. Let v p (x) denote the exponent of the largest 
power of p dividing the integer x and v p (0) = oo. Let S be a set of integers. 
We use (the same) notation v P {S) to denote the expected p- valuation 
of x £ S. If integers in S are random and uniformly distributee^, the 
expected p- valuation v P {S) is 

v P (S) = E [v v {x)\ = Pr(i/„ > 1)+Pr(i/ P > 2) + - • ■ = - + \ + - ■ ■ = — l —. 

xes p p p — 1 

Thus, in an informal (logarithmic) sense, an integer s in S contains an 
expected power p L ^ p ~ v> . 

Let now S be a set of polynomial values f(x). We use (the same) nota- 
tion v p (S) (or v p (f)) to denote the expected p- valuation of the polynomial 
values S. Hensel's lemma gives conditions when a root of / (mod p e ) can 
be lifted to a root of / (mod p e+1 ). 

Lemma 2.2 (Hensel's lemma). Let ri be a root of f(x) modulo an odd 
prime p. 

1. If n is a simple root, f{x) (mod p e ) has an unique root r e = r\ 
(mod p) for each e > 1. 

2. If r e is a multiple rooQ of f (mod p e ) for e > 1, there are two 
possible cases. If p e+1 | f(r e ), then\/i£ [0,p), p e+1 \ f(r e + ip e ). If 
p e+1 \ /( r e), r e cannot be lifted to a root modulo p e+1 . 

Assume now that the integers x leading to the values f{x) £ S are 
uniformly random. There are two cases. First, suppose p \ A, the dis- 
criminant of f(x). p is an unramified prime. Then f(x) (mod p) has only 
simple roots. Let n p be the number of roots. The expected p-valuation of 
polynomial values is u p (f) = n p /(p — 1) (apply the formula above, using 
Pr(fp > e) = n p /p e ). 

The second case is when p | A. Here one may get multiple roots. The 
expected p-valuation may be obtained by counting the number of lifted 
roots. 

In the number field sieve, we want to know the expected p-valuation of 
homogeneous polynomial values F(a,b), where (a, b) is a pair of coprime 
integers, and F(x, y) is the homogenous polynomial corresponding to f(x). 
We assume in the following that (a, b) is a uniformly random pair of 
coprime integers. We have 

v p {F(a,b)) = v p {F{\a,\b)) (3) 

for any integer A coprime to p. A pair of coprime integers (a, b) maps 
to a point (a : 6) on the projective line P (F p ). Because of property Q 
above, pairs for which v p (F(a,b)) > correspond to the points of the 
zero-dimensional variety on P 1 (F P ) defined by the polynomial F. 

1 We consider integer random variables within a large enough bounded sample space. 
2 We say that r e is a multiple root of / (mod p e ) if f'(r e ) = (mod p). 
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The projective line P 1 (F p ) has p+ 1 points, consisting of p affine points 
which can be represented as (a; : 1) with x £ F p , together with the point 
at infinity (1 : 0). Among these, the zeroes of F correspond, for affine 
points (x : 1), to affine roots x £ F p of the dehomogenized polynomial /. 
The point at infinity is a zero of F if and only if the leading coefficient Cd 
of / cancels modulo p. If F has a total of n p affine and projective zeroes 
in P 1 (F P ), then F(a, b) for coprime (a, b) is divisible by p with probability 
n p /(p + l). 

It is also possible to look at (a, b) modulo a prime power p e . Then (a, b) 
maps to an equivalence class (a : b) on the projective line over the ring 
Z/p e Z. The p- valuation of F at (a : b) £ P 1 (Z/p e Z) (an integer between 
and e — 1, or "e or more") conveys the information of what happens 
modulo p e . There are p e + p e_1 points in P 1 (Z/p e Z) (p e affine points of 
the form (x : 1), while the remaining p e ~ 1 points at infinity are written 
as (1 : py)). A coprime pair (a, 6) chosen at random maps therefore to a 
given point in P 1 (Z/p e Z) with probability l/(p e_1 (p + 1)). 

Given an unramified p, let F(x, y) (mod p) have n p affine and pro- 
jective roots (zeroes on P 1 (F P )). In application of the Hensel Lemma 
(applied to / at an affine root x, or to p d ~ 1 f(^) above the possible pro- 
jective root), there is a constant number n p of points (a : b) = P 1 (Z/p e Z) 
such that f p (F(a : b)) > e, as e grows. The expected p-valuation v p (F) is 
thus: 

oo 

For ramified p, simply counting the number n p of affine and projective 
roots modulo p is not sufficient to deduce v p (F). More careful compu- 
tation is needed modulo prime powers, which is addressed in Sections 0] 
and [521 

Murphy [141 p. 49] defines the a(F) function to compare the cumu- 
lative expected p-valuation of polynomial values to random integers of 
similar size. a(F) can be considered as the logarithmic benefit of using 
polynomials values compared to using random integers. 

p<B ' 
V prime 

where the summand rewrites as ( 1 Up ^ ] _2fLE when p is unramified. 

V p+ij p-i 

In the number field sieve, a(F) is often negative since we are interested 
in the case when F(x, y) has more than one root. 

2.4 Steps in polynomial selection 

Polynomial selection can be divided into three steps: polynomial genera- 
tion, size optimization and root optimization. 

In the polynomial selection, we first generate good polynomials in 
terms of the size property. Two efficient algorithms are given by Kleinjung 
[21 [TO]. Once we have generated some polynomial pairs (f(x) = g(x) = 
m,2X — mi) of relatively good size, the size and root properties of these 
polynomials can be further optimized using translation and rotation. 
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• Translation of f(x) by k gives a new polynomial fk{x) denned by 
fk(x) = f(x + k). The root of fk{x) is mi/rri2 — k (mod n). The 
linear polynomial gk(x) is m,2X — mi + km,2 . Translation only affects 
the size property. 

• Rotation by a polynomial X(x) gives a new polynomial f\( x )(x) de- 
fined by fx( x )( x ) = f( x ) + A(x) (m22; — mi). The linear polynomial 
is unchanged ffAfaoC 3 -) — sOe) = ?ri2a: — mi. The root is unchanged. 
A(a;) is often a linear or quadratic polynomial, depending on n and 
the skewness of f(x). Rotation can affect both size and root prop- 
erties. 

Given a polynomial pair, translation and rotation are used to find a 
polynomial of smaller (skewed) norm (cf Equation This is called size 
optimization. 

Many polynomials can have comparable size after size optimization. 
We produce and choose the best polynomials in terms of good a-values. 
This requires that the polynomials have many roots modulo small prime 
and prime powers. This step is referred to as root optimization. 

Given f(x) (or F(x,y)), we can use polynomial rotation to find a 
related polynomial f\( x )(x) (or Fx( x )(x,y)) which has a smaller a but 
similar size. Polynomial rotation may also increase the size of trailing 
coefficients. However, if the skewness of the polynomial is large, the size 
property of the polynomial may not be altered significantly. Hence there 
is some room for rotation if the skewness is large. As an indication of this, 
the skewed L°° norm of /, defined as max; \s l ~ d ^ 2 fi\, remains unchanged 
for f\(x) as long as the trailing coefficients of f\(x) do not dominate. 
This is true for the polynomials generated by the algorithm |f 0| . where 
the skewness for the polynomials is likely to be large. 

We discuss some algorithms for root optimization in the following sec- 
tions. 

3 Root sieve 

We focus on root optimization for quintic and sextic polynomials in this 
chapter. Given a polynomial pair (f,g), we want to find a rotated poly- 
nomial with similar size but better root properties. We consider linear 
rotations defined by f u<v (x) = f(x) + (ux + v)g(x). We want to choose 
(u, v) such that fu,v(x) has a small a- value. 

The straightforward way is to look at individual polynomials fu,v(x) 
for all possible (u, v)'s and compare their a- values. This is time-consuming 
and impractical since the permissible bounds on U, V are often huge. 

Murphy |14l p. 84] describes a sieve-like procedure, namely the root 
sieve, to find polynomials with good root properties. It is a standard 
method to optimize the root property in the final stage of polynomial 
selection. We describe Murphy's root sieve in Algorithm [T] Let B be 
the bound for small primes and U, V be bounds for the linear rotation. 
The root sieve fills an array with estimated a-values. The a-values are 
estimated from p- valuation for small primes p < B. Alternatively, it is suf- 
ficient to calculate the summation of the weighted p- valuation v p (F) logp 
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for the purpose of comparison. The idea of the root sieve is that, when r 
is a root of fu,v(x) (mod p e ), it is also a root of f u +i P <> ,u+jp c (x) (mod p e ). 



Algorithm 1: Murphy's root sieve 

input : a polynomial pair /, g; integers U, V, B; 

output: an array of approximated a- values of dimension U x V; 

1 for p < B, p prime do 

2 for e where p e < B do 

3 for x E [0,p e — 1] do 

4 for u E [0,p e - 1] do 

5 compute v in /(x) + (ita; + v)g(x) = (mod p e ); 

6 update v p {fu+ip" ,v+] P ") by sieving; 



In general, the root sieve does not affect the projective roots signifi- 
cantly. It is sufficient to only consider the affine roots' contribution to the 
a-value. In the end, we identify good slots (those with small a-values) 
in the sieving array. For each slot (polynomial), we can compute a more 
accurate a-value with a large bound B and re-optimize its size. 

We consider the asymptotic complexity of Murphy's root sieve. 



AlogBj 

E y E pv (oiD + ^) 
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logB 



We are interested in small primes and hence Bj log B is small. The sieving 
bounds U, V dominate the running-time 0(UV B / log B). 



4 A faster root sieve 

In the root sieve, we identify the number of roots of "rotated" polynomials 
,fu,v(x) for small primes and prime powers. In most cases, the roots are 
simple, and hence their average p- valuation follows Equation Q. There 
is no need to count the lifted roots for them. We describe a faster root 
sieve based on this idea. 

We use the following facts based on Hensel's lemma. Suppose n is 
a simple root of f(x) (mod p). There exists a unique lifted root r e of 
fix) (mod p e ) for each e > 1. In addition, each lifted root r e is a simple 
root of f(x) (mod p). For convenience, we say r e is a simple root of f(x) 
(mod p e ) if f'(r e ) =£ (mod p). 

Let r e be a simple root of a rotated polynomial f u ,v(x) (mod p e ) for 
e > 1. It is clear that fu+ip c ,v+jp B (x) = fu,v{x) (mod p e ) for integers 
It follows that r e is also a simple root of the rotated polynomials 
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fu+ip", v+jpe (x) (mod p e ). Given a simple root r\ of a polynomial f u ,v(%) 
(mod p), the contribution of the root n to v p {F UtV ) is p/ (p 2 — 1). We can 
update the scortjf] for all rotated polynomials fu+ip,v+jp{x) in a sieve. 

If r e is a multiple root of f(x) (mod p e ) for some e > 1, there are two 
possible cases. If f(r e ) = (mod p e+1 ), then VZ G [0,p), f(r e + lp e ) = 
(mod p e+1 ). There are p lifted roots r e +i satisfying r e +i = (r e + Zp e ) = r e 
(mod p e ), VZ G [0,p). In addition, the lifted roots r e +i are multiple since 
/'(»"e+i) = (mod p). On the other hand, if f(r e ) ^ (mod p e+1 ), r e 
cannot be lifted to a root modulo p e+1 . 

Let r e be a multiple root of a rotated polynomial fu,v{x) (mod p e ) for 
e > 1. It is also a multiple root for all rotated polynomials f u +ip',v+jp e (x) 
(mod p e ). 

Let r be a fixed integer modulo p. We discuss the case when r is a 
multiple root for some rotated polynomial fu,v{x) (mod p). We see that 
f(r) + (ur + v)g(r) = (mod p) and f'(r) + ug(r) + (ur + v)g'(r) 
i ! (modp). Since (ur + v) = —f(r)/g(r) (mod p), we get ug 2 (r) = 
f(r)g'(r) - f'(r)g(r) (modp). 

Therefore, only 1 in p of it's admits a multiple root at r (mod p). 
For the other it's, we can compute v and update the simple contribution 
p/(p 2 — 1) to slots in the sieve array. If r is a multiple root of f u ,v(x) 
(mod p), we have to lift to count the lifted roots. We discuss the details of 
the lifting method in the following sections. For the moment, we describe 
the improved root sieve in Algorithm [5] 



Algorithm 2: A faster root sieve 

input : a polynomial pair /, g; integers U, V, B; 

output: an array of approximated a- values of dimension U x V; 

l for p < B, p prime do 



2 for x G [0,p— 1] do 

3 compute u such that ug 2 (x) = f(x)g'(x) — f'(x)g(x) (mod p); 

4 for u e [0,p — 1] do 

5 compute v such that f{x) + uxg(x) + vg(x) = (mod p); 

6 if u ^ u; 

7 then 

8 update v p {fu+ip,v+jp) in sieving; 

9 else 

10 lift to count multiple roots of fu.v(x) (mod p e ) such that 
(u, v) = (u, v) (mod p), u, v < p e , p e < B and then sieve; 



Let r — x be fixed in Line 2. In Line 3, we compute u such that r is 
a multiple root of fu, v (x) for some v. If u 7^ w, r is a simple root for this 
u, and some v which will be computed in Line 5. If u — u, fu,v(x) admits 
r as a multiple root. We need to lift (up to degree d) to count the roots. 

3 v p (Fu,v) logp, the contribution of the root r\ (mod p) to a(F u ,v)- 



The running-time to do this is about 
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P prime 

The asymptotic running-time has the same magnitude. In practice, how- 
ever, we benefit of not considering the prime powers. For comparison, 
Murphy's root sieve takes about UV X/ P <bO°S B)/Q°EP) operations, while 
Algorithm [2] takes about UVY^ <B 1 operations. Taking B = 200 for in- 
stance. E P < 2 oo(l°g 200)/(Iogp) « 2705 and £ P < 2 oo 1 = 46. 

5 A two-stage method 

We give a two-stage algorithm for the root optimization. The algorithm is 
motivated by previous work by Gower [6], Jason Papadopoulos (personal 
communication), Stahlke and Kleinjung [18], who suggested to consider 
congruence classes modulo small primes. 

If the permissible rotation bounds U, V are large, the root sieve can 
take a long time for each polynomial. This is even more inconvenient 
if there are many polynomials. We describe a faster method for root 
optimization based on the following ideas. 

A polynomial with only a few roots modulo small prime powers p e is 
less likely to have a good a- value. Therefore, rotated polynomials with 
many (comparably) roots modulo small prime powers are first detected. 
A further root sieve for larger prime powers can then be applied. 

In the first stage, we find a (or some) good rotated pair (uo,vq) 
(mod pi 1 ■ ■ ■ p^™ ) such that the polynomial f UQ ,„ (x) has many roots mod- 
ulo (very) small prime powers p^ 1 , • • • , p^" . Let B a be an upper bound 
for p^" . In the second stage, we apply the root sieve in Algorithm [2] to 
the polynomial fu ,v (x) for larger prime powers up to some bound B. 

5.1 Stage 1 

Given f(x), we want to find a rotated polynomial fu ,v (x) which has 
many roots modulo small primes and small prime powers. Let the prime 
powers be p\ x , ■ ■ ■ , p^™. There are several ways to generate fu ,v {x). 

First, we can root-sieve a matrix of pairs (u, v) of size (112=1 pT ) 2 anc ^ 
pick up the best (u, v) pair(s) as (ito,«o). If the matrix is small, there is 
no need to restrict the bound in the root sieve to be B a . We can use the 
larger bound B. If the matrix is large, however, the root sieve might be 
slow. We describe a faster strategy. 

We first find m (or moreQ) individual polynomials fu^v^Piix) (1 < 
i < m) each of which has many roots modulo small p^ . The values Ui 
and Vi are bounded by p^' . We combine them to obtain a polynomial 
fu ,v (x) (mod n^iPi*) using the Chinese Remainder Theorem. The 
polynomial fu ,v (x) (mod p*) has the same number of roots as the in- 
dividual polynomials fu^vi^iix) (mod p') for 1 < k < e^. Hence the 

4 For each pi, we can generate more than one polynomial. In Stage 2 we consider multi-sets 
of combinations. 



combined polynomial is likely to have many roots modulo small prime 
powers pi 1 , •• ■ ,p%?. 

Individual polynomials. To find individual polynomials f Ui ,vi,pi (x) 

that have many roots modulo small prime powers p\ l , we can root-sieve 
a square matrix pi' x p\ % and pick up the good pairs. 

Alternatively, we use a lifting method together with a p 2 -ary tree data 
structure. This seems to be more efficient when p 2e ' is large. For each 
p\ T , we construct a tree of height d and record good (u, v) pairs during 
the lift. The lift is based on Hensel's lemma. For convenience, we fix f(x) 
(mod p) where p — Pi and e = for some i. We describe the method. 

We create a root node. In the base case, we search for polynomials 
fu,v(x) (mod p) [u, v £ [0, p)) which have many roots and record them 
in the tree. There can be at most p 2 level-1 leaves for the root node. In 
practice, one can discard those leaves with fewer (comparably) roots and 
only keep the best branch. 

Let a level-1 leaf be (u, v) (mod p). A simple root is uniquely lifted. If 
the polynomial f u<v (x) (mod p) only gives rise to simple roots, we already 
know the exact p- valuation of f u , v (x). In case of multiple roots, we need 
to lift and record the lifted pairs. Assume that f u ,v(x) (mod p) has some 
multiple root r m and some simple root r s . We want to update the p- 
valuation for rotated polynomials 

/ e-1 e-1 \ 

/(*)+ (« + E ik P k ) x + { v + E 3kP k ) 9(x) (mod p e ) (5) 

\ fe=l k=l I 

where each ik,jk £ [0,p). We give the following procedure for the lifting. 

1. For a simple root r„, we find out which of the rotated polynomials 
fu+i P .v+ip{x) (mod p 2 ) admit r 3 as a root. If f u +ip,v+jp(r s ) = 
(mod p 2 ) for some then 

(ir s + j)g(r s ) + f u ,v{r s )/p = (mod p). (6) 

Hence the set of (i, j)'s satisfies a linear congruence equation. For 
simple roots, there is no need to compute the lifted root. It is suf- 
ficient to update the p-valuation contributed by r s to polynomials 

fu + ip,v+jp{x) . 

2. Let r m be a multiple root of fu,v(x) (mod p). If a rotated polynomial 
fu+ip,v+j P (x) (mod p 2 ) admits r m as a root for some all the 
{r m + lp} (0 < I < p) are also roots for the polynomial. In addition, 
fL+ip,v+jp{ r m + lp) =0 (modp). We record the multiple roots 
(r m + lp} together with the {(u + ip, v + jp)} pairs. The procedure 
also works for the lift from p e to p e+1 for higher e's. 

We consider the memory usage of the p 2 -ary tree. If r is a root of 
fu,v(x) (mod p), Equation ((6]) shows a node (u, v) (mod p) gives p lifted 
nodes (it + ip,v + jp) (modp 2 ) for some (i,j)'s. Since fu,v(x) (mod p) 
can potentially have other roots besides r, there could be at most p 2 pairs 
(u + ip, v + jp) (mod p 2 ). The procedure also needs to record the multiple 
roots for each node. We are mainly interested in the bottom level leaves 
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of the tree, those (u, v) (mod p e ). It is safe to delete the tree path which 
will not be used anymore. Hence a depth-first lifting method can be used. 
In practice, the memory usage is often smaller than a sieve array of size 

p 2e - 

For each p, we find a polynomial that either has many simple roots or 
many multiple roots which can be lifted further. Tiny primes p's are more 
likely to be ramified. Hence we are more likely to meet multiple roots for 
tiny p. 

CRTs. For each p, we have generated some polynomial(s) rotated by 
(ux+v)g(x) (mod p e ) which have comparably good expected p-valuation. 
For convenience, we identify the rotated polynomial by pair (u,v). 

Stage 1 repeats for prime powers p^ 1 , ■ ■ ■ , p^™ . Let M = YlTLi Pi' ■ We 
generate the multi-sets combinations of pairs {{u, v)} and recover a set 
of {(uo,vo)} (mod M). We fix such a pair (rotated polynomial) (uo,vo) 
(mod M). 

The whole search space is an integral lattice of Z 2 . In Stage 2, we 
want to root-sieve on the sublattice points defined by (wo + 7-M, vo + fiM) 
where (7,/?) £ Z 2 . The sublattice points are expected to give rotated 
polynomials with good root properties, since the polynomials have many 
roots modulo p^ 1 , ■ ■ ■ , p e ^ . 

We often choose the pi's to be the smallest consecutive primes since 
they are likely to contribute most to the a- value. The exponents ei in 
prime powers p\ l need some more inspection. If a is too small, the sieving 
range (7, /3) G 7? can be large. If e* is too large, M is large and hence 
some polynomials which have good size property might be omitted in the 
root sieve. One heuristic is to choose pi' < p? for i > j, i,j £ [1, m]. To 
determine m, one can choose M to be comparable to the sieving bound U. 
Assume that M ~ U. We can discard those (uo,vo)'s such that uo > U. 
If uo is comparable to U, it is sufficient to use a line sieve for constant 
rotations. 

Remark 5.1. In the implementation, we may want to tune the parameters 
by trying several sets of parameters such as various pi's and ej's. We 
can run a test root sieve in short intervals. The set of parameters which 
generates the best score is then used. 

5.2 Stage 2 

In Stage 2, we apply the root sieve in Algorithm[2]to polynomial f uo ,v {x), 
perhaps with some larger prime bound. In the root sieve, one can reuse 
the code from Stage 1, where the updates of a-values can be batched. We 
describe the method as follows. 

Sieve on sublattice. Let M = YliLi Pi* an< i ( u o,"o) be fixed from 
Stage 1. In the second stage, we do the root sieve for (larger) prime powers 
on the sublattice defined by {(uo + jM, vo + /3M)} where 7, j3 £ Z. Let p 
be a prime and (mod p k ) be a root of 

f(x) + ({uo + jM)x + {vo + pM)^jg{x) (mod p k ) 
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for some fixed integers 7, /3. The sieve on the sublattice follows from 

/(V fc )+((«o + M( 7 + ip fc ))r fc + (w + M(j8 + jp fe ))) <?(r fe ) = (mod p fe ) 

for integers i,j £ Z. We consider the root sieve for a fixed prime p in 
Algorithm 

Let f(x),g(x),M,uo,vo be fixed from Stage 1. In Algorithm [2] we 
assume u, r are fixed for the moment. Let p be a prime not dividing M. 
The sieve array has approximate size \U/M\ x \V/M\. Each element 
(7,/?) in the sieve array stands for a point (uo + jM,vo + /3M) in Z 2 . 
We solve for v in f(r) + urg(r) + vg(r) = (mod p). Knowing (u, v), we 
can solve for (7, 0) in u = uo + yM (mod p) and v = Vo + /3M (mod p), 
provided that p\M. 

For the moment, we fix integers 7, j3. If r is a simple root, it is sufficient 
to sieve (7 + ip, {3 + jp) for various (i, j)'s and update the p- valuation 
p log p/(p 2 — 1) to each slot. If r is a multiple root, we can use a similar 
lifting procedure as in Stage 1. We describe the recursion to deal with 
multiple roots in Algorithm [21 



Algorithm 3: Recursion for multiple roots 

input : a polynomial pair /, g; integers U, V, B; node (u,v), tree height 

e, current level k, prime p; 
output: updated a-values array; 

1 for multiple roots r of f UjV (x) (mod p) do 

2 for k < e do 

3 compute (i,j)'s in (ir + j)g(r) + fu,v(f)/p = (mod p); 

4 create child nodes (u + ip k , v + jp k ) with roots 
{r + lp k }, VI e [0,p); 

5 recursively call Algorithm [3] on (u,v)'s leftmost child node; 

6 change coordinates for current node (u, v) and sieve; 

7 delete current node and move to its sibling node or parent node; 



From Stage 1, we know uo,vo- In Algorithm [2l we fix u,r and solve 
for v. Given a multiple root r of f(x) (mod p*), we find pairs (u',v') 
such that f u '. v '(r) = (mod p k+1 ) where u' = u (mod p k ) and v' = v 
(mod p k ). We can construct nodes representing the («', v') pairs together 
with their roots. In the recursion, we compute the lifted nodes in a depth- 
first manner. Once the maximum level p e is reached, we do the root sieve 
for the current nodes and delete the nodes which have been sieved. 

When a lifted tree node (u',v') (mod p fe ) is created, the number of 
roots for f u > !V >(x) (mod p k ) is known. In the root sieve, the Q-scores can 
be updated in a batch for all the roots of f u '. v '(x) (mod p k ). For each 
node (u',v') : we also need to compute the corresponding coordinates in 
the sieve array. 

Primes p dividing M . We have assumed that p is a prime not divid- 
ing M. From Stage 1, M is a product of prime powers p?* for 1 < i < m. 
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For accuracy, we can also consider primes powers p^ with e' t 7^ ej such 
that pi appears in the M. Let r be root of fu,v{x) (mod p). If r is a 
simple root, there is no need to consider any liftings. Hence we consider 
polynomials fu,v{x) (mod p) which have a multiple root. 

We fix some p — pi and e = et, which are used in Stage 1. Let u,v,p 
be fixed in Algorithm [2] Let e' be the exponent of p that we want to 
consider in Stage 2. There are two cases depending on e'. 

If e' < e, the points on the sublattice have equal scores contributed by 
roots modulo p e . It is sufficient to look at the multiple roots modulo p k 
for k < e . In Algorithm [2] we either sieve all slots of the array or do not 
sieve at all. Given u, v,p, k, if v = «o (mod p k ) in v = v + /3M (mod p k ), 
we need to sieve the whole array. This can be omitted because it will 
give the same result for each polynomial and we only want to compare 
polynomials. If vo ^ v (mod p k ), no slot satisfies the equation. Therefore, 
it is safe to skip the current iteration when e' < e. 

If e' > e, the rotated polynomials (u, v) (mod p k ) for e < k < e' may 
have different behaviors. We describe some modifications in the lifting 
procedure. Let u = uo (mod p fe ), r, vo be fixed in Algorithm [2] We 
compute v. We want to know which points (polynomials) on the sieve 
array are equivalent to (u, v) (mod p k ). 

For k < e, the situation is similar to the case when e' < e. If the 
equation v = vo (mod p k ) is satisfied, we record the node (u, v) (mod p k ) 
for further liftings. There is no need to sieve since all slots on the sieve 
array have equal scores for roots modulo p k . If Vo ^ v (mod p) where 
k = 1, we have neither to root-sieve nor record the node. Let e < k < e . 
If (v! , v') (mod p k ) satisfies u' = u (mod p) and v' = v (mod p), we need 

vq + fiM = v' (mod p ). 

The equation is solvable for f3 only if 

vq = v' (mod p k ). 

Hence it is safe to discard those (it, v) (mod p) such that u = uo (mod p) 
but 

v ^ vo (mod p) . 

On the other hand, we consider some k in e < k < e' . In the lifting 
procedure, we record nodes without sieving until we reach the level-(e + 1) 
nodes. Starting from a node (u, v) modulo p e+1 , that is k > e, we want 
to solve the equation 

wo + PM = v (mod p k ). 
The depth-first lifting procedure shows that 

vo = v (mod p e ). 
Hence /3 is solvable in the following equation 

vq-v M fc _ e 
h — P = (mod p ) 

since gcd(M/p e ,p) = 1. In the root sieve, we step the array by /3 + jp k ~ e 
for various j. 
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5.3 Further remarks and improvements 



Let (U, V) be the rotation bounds for the polynomial. The root sieve in 
Algorithm [2] runs asymptotically in time UVB/\ogB (ignoring constant 
factors). In Stage 2, the searching space is restricted to a sublattice de- 
termined by M = n^LiPi'i where the parameters pi's depend on Stage 
1. Hence, the root sieve in Stage 2 runs in time about UVB/(M 2 log B). 

In Stage 2, the points not on the sublattice are discarded since com- 
pared to points on the sublattice they have worse p- valuation for those p's 
in Stage 1. We assumed that they were unlikely to give rise to polynomials 
with good root properties. However, a polynomial could have good a(F) 
while some p in M gives a poor p-valuation. This often happens when 
some p'-valuation of p' \ M, those ignored in Stage 1, is exceptionally 
good, and hence mitigates some poor p-valuation where p | M. 

Alternatively, we can use a root sieve to identify good rotations in 
Stage 1 for some small sieving bounds (U',V'). Then we examine the 
pattern of p-valuation of these polynomials and decide the congruence 
classes used in Stage 2. 

We have ignored the size property of polynomials in the algorithms. 
We have assumed that polynomials rotated by similar (u, v)'s have com- 
parable size. In practice, some trials are often needed to decide the sieving 
bounds (U, V). We give some further remarks regarding the implementa- 
tion. 

Block sieving. The root sieve makes frequent memory references to 
the array. However, there is only one arithmetic operation for each array 
element. The time spent on retrieving memory often dominates. For 
instance, the root sieve may cause cache misses if the sieve on p steps over 
a large sieve array. A common way to deal with cache misses is to sieve 
in blocks. 

We partition the sieving region into multiple blocks each of whose size 
is at most the cache size. In the root sieve, we attempt to keep each block 
in the cache while many arithmetic operations are applied. The fragment 
of the block sieving is described in Algorithm [4] 



Algorithm 4: Block sieving 

input : a polynomial pair /, g; integers U, V, B; 

output: an array of approximated a- values in dimension U x V; 

1 for x < B do 

2 for each block do 

3 for p where x < p < B do 

4 



We have also changed the order of iterations to better facilitate the 
block sieving. This might give some benefits due to the following heuristic. 
In Algorithm [2] when p is small, polynomial roots x modulo p are small. 
The number of roots x < p blocked for sieving is also limited. Instead 
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we block primes p. If x is small, there are still many p's which can be 
blocked. 

For multiple roots, we might need to sieve in steps p k for k > 1. When 
p k is not too small, each block has only a few (or none) references. In 
this case, we may use a sorting-based sieving procedure like the bucket 
sieve Q]. 

Arithmetic. The coefficients of the rotated polynomials are multiple 
precision numbers. Since p e can often fit into a single precision integer, it 
is sufficient to use single precision in most parts of the algorithms. 

The algorithms involve arithmetic on p k for all k < e. It is sufficient 
to store polynomial coefficients modulo p e and do the modulo reduction 
for arithmetic modulo p k . Let D be a multiple precision integer. In the 
algorithm, we use a single precision integer S instead of D where S = D 
(mod p e ). If x = D (mod p k ) for k < e, it is clear that x = S (mod p k ). 
Hence we can use the S in the root optimization. 

In addition, the range of possible a- values is small. We may use short 
integers to approximate the a- values instead of storing floating point num- 
bers. This might save some memory. 

Quadratic rotation. Sextic polynomials have been used in the factor- 
izations of many large integers such as RSA-768. Rotations by quadratic 
polynomials can be used for sextic polynomials if the coefficients and skew- 
ness of the polynomials are large. We have assumed that W is small in 
fw,u,v{x) and we restricted to use linear rotations in this section. If the 
permissible bound for W is large, we can use a similar idea to that in 
Stage 1 to find good sublattices in three variables. At the end of Stage 1, 
a set of polynomials having good a-values are found which are defined 
by rotations of (wo, uo, vo)'s. In Stage 2, we root-sieve on the sublattice 
{{w + SM, u + jM, v + I3M)} where 8, 7, /9 € Z. 

6 Conclusion 

Root optimization aims to produce polynomials that have many roots 
modulo small primes and prime powers. We gave some faster methods 
for root optimization based on Hensel's lifting lemma and root sieve on 
congruences classes modulo small prime powers. The algorithms described 
here have been implemented and tested in practice. The implementation 
can be found in C ADO-NFS [2]. 
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